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1) ■ Abstract 

In this paper we address the issue of modeling electricity loads. After analyzing 
properties of the deseasonalized loads from the California power market we fit an 
& ARMA(1,6) model to the data. The obtained residuals seem to be independent 

but with tails heavier than Gaussian. It turns out that the hyperbolic distribution 
provides an excellent fit. 
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1 Introduction 



During the last decade we have witnessed radical changes in the structure of electricity 
markets world-wide. For many years it was argued convincingly that the electricity in- 
dustry was a natural monopoly and that strong vertical integration was an obvious and 
efficient model for the power sector. However, recently it has been recognized that compe- 
tition in generation services and separation of it from transmission and distribution would 
be the optimal long-term solution. Restructuring has been designed to foster competition 
and create incentives for efficient investment in generation assets [7,12,16]. 

While the global restructuring process has achieved some significant successes, serious 
problems - some predictable, others not - have also arisen. The difficulties that have 
appeared were partly due to the flaws in regulation and partly to the complexity of the 
market. 

When dealing with the power market we have to bear in mind that electricity cannot 
simply be manufactured, transported and delivered at the press of a button. Moreover, 
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electricity is non-storable, which causes demand and supply to be balanced on a knife- 
edge. Relatively small changes in load or generation can cause large changes in price and 
all in a matter of hours, if not minutes. In this respect there is no other market like it. 

Californians are very well aware of this. In January 2001 California's energy market was 
on the verge of collapse. Wholesale electricity prices have soared since summer 2000, see 
the top panel of Fig. 1. The state's largest utilities were threatening that they would 
be bankrupted unless they were allowed to raise consumer electricity rates by 30%; the 
California Power Exchange suspended trading and filed for Chapter 11 protection with the 
U.S. Bankruptcy Court. How could this have happened when deregulation was supposed 
to increase efficiency and bring down electricity prices? It turns out that the difficulties 
that have appeared are intrinsic to the design of the market, in which demand exhibits 
virtually no price responsiveness and supply faces strict production constraints [4,13]. 

Another flaw of deregulation was the underestimation of the rising consumption of elec- 
tricity in California. The soaring prices and San Francisco blackouts clearly showed that 
there is a need for sophisticated tools for the analysis of market structures and modeling 
of electricity load dynamics [2,9]. In this paper we investigate whether electricity loads in 
the California power market can be modeled by ARMA models. 



2 Preparation of the data 



The analyzed database was provided by the University of California Energy Institute 
(UCEI, www.ucei.org). Among other data it contains system- wide loads supplied by Cal- 
ifornia's Independent (Transmission) System Operator. This is a time series containing 
the load for every hour of the period April 1st, 1998 - December 31st, 2000. Due to a very 
strong daily cycle we have created a 1006 days long sequence of daily loads. Apart from 
the daily cycle, the time series exhibits weekly and annual seasonality, see the bottom 
panel of Fig. 1. Because common trend and seasonality removal techniques do not work 
well when the time series is only a few (and not complete, in our case ca. 2.8 annual 
cycles) cycles long, we restricted the analysis only to two full years of data, i.e. to the 
period January 1st, 1999 - December 31st, 2000, and applied a new seasonality reduction 
technique [15]. 

The seasonality can be easily observed in the frequency domain by plotting a sample 
analogue of the spectral density, i.e. the periodogram 



I n {uk) = ~ 

n 



J2x t exp{-2ni(t - l)u k } 
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where {xi, ...,x n } is the vector of observations, uj^ = k/n, k — 1, ..., [n/2] and [x] denotes 
the largest integer less then or equal to x. In the top panel of Fig. 2 we plotted the 
periodogram for the system-wide load. It shows well-defined peaks at frequencies corre- 
sponding to cycles with period 7 and 365 days. The smaller peaks close to uuk = 0.3 and 0.4 
indicate periods of 3.5 and 2.33 days, respectively. Both peaks are the so called harmonics 
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Fig. 1. California Power Exchange daily average market clearing prices (top panel) and California 
power market daily system-wide load (bottom panel) since January 1st, 1999 until December 31st, 
2000. The annual and weekly seasonalities are clearly visible. 

(multiples of the 7-day period frequency) and indicate that the data exhibits a 7-day pe- 
riod but is not sinusoidal. The weekly period was also observed in lagged autocorrelation 
plots [14]. 

To remove the weekly cycle we used the moving average technique [5]. For the vector of 
daily loads {xi, ..., #731} the trend was first estimated by applying a moving average filter 
specially chosen to eliminate the weekly component and to dampen the noise: 



rh t = -(x t -3 + ■■■ + ^+3), 



(2) 



where t = 4, ...,728. Next, we estimated the seasonal component. For each k = 1,...,7, 
the average Wk of the deviations {(xk+7j — rhk+7j),4: < k + 7j < 728} was computed. 
Since these average deviations do not necessarily sum to zero, we estimated the seasonal 
component Sk as 
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Sk = w k 
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where k — 1, ..., 7 and §k = Sk-7 for k > 7. The deseasonalized (with respect to the 7-day 
cycle) data was then defined as 



d t = x t — s t for t — 1, ..., 731. 



(4) 
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Fig. 2. Periodogram of the California power market daily system-wide load since January 1st, 
1999 until December 31st, 2000 {top panel). The annual and weekly frequencies are clearly visible. 
Periodogram of the load returns after removal of the weekly and annual cycles {bottom panel). 
No dominating frequency can be observed. 

Finally we removed the trend from the deseasonalized data {tit} by taking logarithmic 
returns r t = log(d t +i/d t ), t — 1, ..., 730. 

After removing the weekly seasonality we were left with the annual cycle. Unfortunately, 
because of the short length of the time series (only two years), the method applied to 
the 7-day cycle could not be used to remove the annual seasonality. To overcome this we 
applied a new method which consists of the following [15]: 

(i) calculate a 25-day rolling volatility [8] for the whole vector {r\, ...,r73o}; 
(ii) calculate the average volatility for one year, i.e. in our case 



v t 



^1999 I ^2000 



(5) 



(iii) smooth the volatility by taking a 25-day moving average of vf, 

(iv) finally, rescale the returns by dividing them by the smoothed annual volatility. 

The obtained time series (see the top panel of Fig. 3) showed no apparent trend and 
seasonality (see the bottom panel of Fig. 2). Therefore we treated it as a realization 
of a stationary process. Moreover, the dependence structure exhibited only short-range 
correlations. Both, the autocorrelation function (ACF) and the partial autocorrelation 
function (PACF) rapidly tend to zero (see the bottom panels of Fig. 3), which suggests 
that the deseasonalized load returns can be modeled by an ARMA-type process. 
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Fig. 3. Load returns after removal of the weekly and annual cycles (top panel). The ACF (bottom 
left panel) and PACF (bottom right panel) for the mean-corrected deseasonalized load returns. 
Dashed lines represent the bounds ±1.96/\/730, i.e. the 95% confidence intervals of Gaussian 
white noise. 

3 Modeling with ARMA processes 



The mean-corrected (i.e. after removing the sample mean=0. 0010658) deseasonalized load 
returns were modeled by ARMA (Autoregressive Moving Average) processes 



X t — <f>iX t -i — ... — 4>pX t ^ p — Z t + Q\Z t -\ + ... + Q q Zt- q , t — 1, ...,n, 



(6) 



where (p, q) denote the order of the model and {Z t } is a sequence of independent, iden- 
tically distributed variables with mean and variance a 2 (denoted by iid(0, a 2 ) in the 
text). 



The maximum likelihood estimators 



■ (0i, ..., <f)p), 9 — ($i, ..., 9 q ) and a 2 of the parame- 
ters = (0i, ..., P ), 6 = ($i, ..., 9 q ) and a 2 , respectively, were obtained after a preliminary 
estimation via the Hannan-Rissanen method [5] using all 730 deseasonalized returns. The 
parameter estimates and the model size (p, q) were selected to be those that minimize 
the bias-corrected version of the Akaike criterion, i.e. the AICC statistics 



AICC 



-21nL + 



2(p + q+l)n 

n — p — q — T 



(7) 



where L denotes the maximum likelihood function and n = 730. 
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Fig. 4. The residuals obtained from the ARMA(1,6) model (top panel). The ACF (bottom left 
panel) and the PACF (bottom right panel) of the residuals. Dashed lines represent the bounds 
±1.96/^/730. 

The optimization procedure led us to the following ARMA(1,6) model (with 64 = #5 = 0) 



X t = 0.332776X t _i + 

+Z t - 0.383245Z t _i - 0.12908Z t _ 2 - 0.149307Z t _ 3 - 0.0531862Z t _ 6 , (8) 

where t — 1, ..., 730 and {Z t } ~ iid(0, 0.838716). The value of the AICC criterion obtained 
for this model was AICC=1956.294. 

In order to check the goodness of fit of the model to the set of data we compared the 
observed values with the corresponding predicted values obtained from the fitted model. 
If the fitted model was appropriate, then the residuals 



W t 



x.-x, 



HVV 



t= 1,...,730, 



(9) 



ft-il 



where X t (<f>, 9) denotes the predicted value of X t based on Xi, ..., X t _i and q t _i = E(X t — 
X t ) 2 /a 2 , should behave in a manner that is consistent with the model. In our case this 
means that the properties of the residuals should reflect those of an iid noise sequence 
with mean and variance a 2 . 

The residuals obtained from the ARMA(1,6) model fitted to the mean-corrected deseason- 
alized load returns are displayed in the top panel of Fig. 4. The graph gives no indication 
of a nonzero mean or nonconstant variance. The sample ACF and PACF of the residuals 



Table 1 

Test statistics and p-values for the residuals. 

Test Test statistics value p-value 

Portmanteau 15.03 (0.7747) 

Turning point 464 (0.0609) 

Difference-sign 361 (0.6536) 

Rank 131090 (0.5529) 



fall between the bounds ±1.96/v730 indicating that there is no correlation in the series, 
see the bottom panels of Fig. 4. Recall that for large sample size n the sample autocor- 
relations of an iid sequence with finite variance are approximately iid with distribution 
N(Q,l/n). Therefore there is no reason to reject the fitted model on the basis of the 
autocorrelation or partial autocorrelation function. However, we should not rely only on 
simple visual inspection techniques. For our results to be more statistically sound we per- 
formed several standard tests for randomness. The results of the portmanteau, turning 
point, difference-sign and rank tests are presented in Table 1. Short descriptions of all 
applied tests can be found in the Appendix. 

As we can see from Table 1, if we carry out the tests at commonly used 5% level, the tests 
do not detect any deviation from the iid behavior. Thus there is not sufficient evidence to 
reject the iid hypothesis. Moreover, the order p = of the minimum AICC autoregressive 
model for the residuals also suggests the compatibility of the residuals with white noise, 
see the Appendix. Therefore we may conclude that the ARMA(1,6) model (defined by eq. 
(8)) fits the mean-corrected deseasonalized load returns very well. 



4 Distribution of the residuals 



In the previous Section we showed that the residuals are a realization of an iid(0,a 2 ) se- 
quence. But what precisely is their distribution? The answer to this question is important, 
because if the noise distribution is known then stronger conclusions can be drawn when 
a model is fitted to the data. There are simple visual inspection techniques that enable 
us to check whether it is reasonable to assume that observations from an iid sequence are 
also Gaussian. The most widely used is the so-called normal probability plot, see Fig. 5. 
If the residuals were Gaussian then they would form a straight line. Obviously they are 
not Gaussian - the deviation from the line is apparent. This deviation suggests that the 
residuals have heavier tails. 

However, we have to bear in mind that in order to comply with the ARMA model as- 
sumptions the distribution of the residuals must have a finite second moment. In the class 
of heavy-tailed laws with finite variance the hyperbolic distribution seems to be a natural 
candidate. 

The hyperbolic law was introduced by Barndorff-Nielsen [1] for modeling the grain size 
distribution of windblown sand. It was also found to provide an excellent fit to the distri- 
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Fig. 5. The normal probability plot of the residuals obtained from the ARMA(1,6) model. If the 
residuals were Gaussian then they would form a straight line. 

butions of daily returns of stocks from a number of leading German enterprises [6,10]. The 
name of the distribution is derived from the fact that its log-density forms a hyperbola. 
Recall that the log-density of the normal distribution is a parabola. Hence the hyperbolic 
distribution provides the possibility of modeling heavy tails. 

The hyperbolic distribution is denned as a normal variance-mean mixture where the 
mixing distribution is the Inverse Gaussian law. More precisely, a random variable has 
the hyperbolic distribution if its density is of the form 



f{x;a,(3,5,fi) 



V^^F 



2a5K 1 (5y/a 2 - ft 2 



exp\-aJ5 2 + (x-i^) 2 + (3(x-n)\, (10) 



where the normalizing constant Ki(i) = | / °° exp{ — \t{x + -)}dx, t > 0, is the modified 
Bessel function with index 1, the scale parameter S > 0, the location parameter \x G R 
and < \/3 \ < a. The latter two parameters - a and j3 - determine the shape, with a 
being responsible for the steepness and (3 for the skewness. 

Given a sample of independent observations all four parameters can be estimated by the 
maximum likelihood method. In our studies we used the 'hyp' program [3] to obtain the 
following estimates 



a 



1.671304, = -0.098790, 8 = 0.298285, fi = 0.076975. 



The empirical probability density function (PDF) - to be more precise: a kernel estimator 
of the density - together with the estimated hyperbolic PDF are presented in Fig. 6. 
We can clearly see that, on the semi- logarithmic scale, the tails of the residuals' density 
form straight lines, which justifies our choice of the theoretical distribution. The adjusted 
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Fig. 6. The empirical probability density function (a kernel estimator of the density) and the 
approximating hyperbolic PDF on the semi-logarithmic scale. 

Kolmogorov statistics K = T/nsup x \F(x) — F n (x)\, where F(x) is the theoretical and 
F n (x) is the empirical cummulative distribution function, returns the value K = 1.5652. 
This indicates that there is not sufficient evidence to reject the hypothesis of the hyberbolic 
distribution of the residuals at the 1% level. For comparison we fitted a Gaussian law to 
the residuals as well. In this case the adjusted Kolmogorov statistics returned K = 1.8019 
causing us to reject the Gaussian hypothesis of the residuals at the same level. 



5 Conclusions 



Due to limited monitoring in a power distribution system its loads usually are not known 
in advance and can only be forecasted based on the available information. In this paper 
we showed that it is possible to model deseasonalized loads via ARMA processes with 
heavy-tailed hyperbolic noise. This method could be used to forecast loads in a power 
market. Its effectiveness, however, still has to be tested and will be the subject of our 
further research. 



Appendix: Tests for randomness 



The portmanteau test. Instead of checking to see if each sample autocorrelation p(j) 
falls inside the bounds ±1.96/^/n, where n is the sample size, it is possible to consider a 
single statistic introduced by Ljung and Box [11] Q = n{n + 2) Y%=i f> 2 (j)/( n ~j)i whose 
distribution can be approximated by the \ 2 distribution with h degrees of freedom. A 



large value of Q suggests that the sample autocorrelations of the observations are too 
large for the data to be a sample from an iid sequence. Therefore we reject the iid 
hypothesis at level a if Q > Xi- a W, where Xi- a * s the (1 — a ) quantile of the x 2 
distribution with h degrees of freedom. 

The turning point test. If y 1; ...,y n is a sequence of observations, we say that there is 
a turning point at time i (1 < i < n) if yi-i < yi and yi > yi+\ or if y^\ > yi and 
yi < yi+i- In order to carry out a test of the iid hypothesis (for large n) we denote the 
number of turning points by T (T is approximately N(fj,x, 0r), where \xt = 2(n — 2)/3 
and erf, = (16n— 29)/90) and we reject this hypothesis at level a if |T— Ht\/&t > ^i-a/2, 
where $i_ Q /2 is the (1 — a/2) quantile of the standard normal distribution. The large 
value of T — [it indicates that the series is fluctuating more rapidly than expected for an 
iid sequence; a value of T — \xt much smaller than zero indicates a positive correlation 
between neighboring observations. 

The difference-sign test. For this test we count the number S of values % such that 
yi > yi-i, i = 2, ..., n. For an iid sequence and for large n, S is approximately N(/j,g, a s ), 
where jig = (n — 1)/2 and a s = (n+l)/12. A large positive (or negative) value of S — /j,s 
indicates the presence of an increasing (or decreasing) trend in the data. We therefore 
reject the assumption of no trend in the data if \S — [As\/&s > ®i-a/2- 

The rank test. The rank test is particularly useful for detecting a linear trend in the 
data. We define P as the number of pairs (i,j) such that y,j > y; L and j > i, i — 
1, ...,n — 1. For an iid sequence and for large n, P is approximately N(fip,<7p), where 
Up = n(n — l)/4 and up = n(n — l)(2n + 5)/72. A large positive (negative) value 
of P — up indicates the presence of an increasing (decreasing) trend in data. The iid 
hypothesis is therefore rejected at level a if \P — fj,p\/ap > $1-0/2- 

The minimum AICC AR model test. A simple test for whiteness of a time series is 
to fit autoregressive models of orders p = 0, 1, ...,p max , for some large p ma x, an d to 
record the value of p for which the AICC value attains the minimum. Compatibility of 
these observations with white noise is indicated by selection of the value p — 0. 
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